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SUMMARY 

A summary of fluid/structure interaction capabilities for the NASTRAN 
computer program is presented. The paper concentrates on indirect applications 
of the program towards solving this class of problem; for completeness and 
comparitive purposes, direct usage of NASTRAN will be briefly discussed. The 
solution technology addresses both steady state and transient dynamic response 
problems. 


INTRODUCTION 

A substantial amount of activity is in progress in the general area of 
applying the NASTRAN computer program to fluid /structure interaction problems. 
The class of problems under consideration is limited to linear elastic struc- 
tures in contact with a fluid. The fluid constitutive equation is represented 
by an acoustical type small deformation relation wherein pressure is propor- 
tional to the divergence of the displacement field (per cent change in elemental 
volume) . The time domain character, of the problems treated are either transient 
(usually incident step pressure waves with an exponential decay wave form) or 
steady state (acoustical induced response resulting from a harmonic train of 
incident or radiating pressure waves) . 

The direct use of NASTRAN to solve problems in the category described 
above is documented in the NASTRAN program manuals, therefore the paper will 
only briefly mention direct usage for completeness and comparitive purposes. 
Instead, the paper concentrates on nonstandard flu id /structure applications of 
NASTRAN that range anywhere from employing the program directly (through an 
analogy argument) to using the program capability indirectly (in conjunction 
with auxiliary post-processing programs) . 

Special attention is given to the case where an elastic structure is 
completely submerged in a limitless fluid domain. Five methods are presented 
for handling the modeling problem of having to represent an infinite fluid 
region with only a finite number of elements. Methodology is covered that 
enables one to either eliminate the need for any fluid elements at all (through 
the proper handling of the fluid /structure interface) or requires one to only 
model a tractable finite number of fluid elements. Two of the five methods 
cover transient problems and the remaining ones are for steady state problems. 
Exact versus finite element solutions are presented for most of the methods 
covered in the paper. 
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STEADY STATE-UNCOUPLED FLUID FIELD 


This class of problems treats the case where a totally submerged elastic 
structure is interacting with an infinite fluid domain. The structure is 
represented by finite elements whereas fluid field is represented by a continuum. 
This category of problem is not handled by any of the rigid formats currently 
in NASTRAN. A special post-processing program called FIST (reference 4) was 
written that will accept basic mode shape information directly from the NASTRAN 
rigid format output tape. FIST is designed to process this information into 
the desired solution for the complete fluid/structure interaction response. 

The method centers about the process of obtaining a relation between the 
fluid /structure interface fluid pressure and the interface normal velocity. 

Once this relationship is determined, the uncoupling process unfolds. The 
starting point for formulations of this type (references 6-9) is the Helmholtz 
integral, ref erence 5, where for any point, x, on the closed submerged surface, 

S, which interfaces with the fluid, the total pressure p(x) on the surface is 
related to the normal velocity, w, on the surface by the integral relation 

p(x) = p 1 (x) - lj p(y) ~ ds(y) + UbipJ" w(y)G(x,y)dS(y) (4) 


where y is a dummy variable for any position x e S, p the mass density of the 
fluid and G is the free space Green’s function given by 


G(x,y) 


exp (-ioj 1 x-y j /c) 
4 tt x-y 


(5) 


The development to follow in this subsection on harmonic analysis follows 
reference 2 (modified for incident pressure by the method given in reference 3) 
for the first part of this subsection on a direct solution to the problem and 
follows reference 4 for the modal solution to the problem. The partial deriv- 
ative of G with respect to n(y)_denotes_the rate of change of G in the direction 
normal_to the surface at point y, and |x-y| denotes the distance between the 
x and y points. 


The next step is to obtain a discrete version of equation (4) which is 
accomplished by representing the surface pressure and normal velocity in terms 
of a linear combination of scalar basis functions ^ defined as 


p(x) 



P n *„M 


( 6 ) 


w(x) 


N 


where N denotes the number of surface grid points in contact with the fluid. 
For example, reference 4 has used a cubic spacial distribution con- 
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sis tent with the finite element displacement fields for the neighboring struc- 
tural elements. This is in contrast to reference 10 which employs a piecewise 
constant distribution of pressure over the interface zones of the structure or 
to reference 11 which employs a quadratic distribution. Employing a higher 
order basis function has the advantage that the same solution accuracy can be 
achieved with a coarser interface mesh. This fact ultimately results in a cost- 
effective computer program that should run more efficiently, while maintaining 
the same accuracy, when employing the higher order distribution basis functions. 

Upon substituting equations (6) into equation (4) and evaluating equation 
4 over a discrete set of points (xj , j = 1, 2, J) corresponding to the fluid 
structure interface node points, one obtains 

[L] {P> = [rHw 1 } + {P 1 } (7) 

where [L] and [R] are JxJ matrices and {P 1 } is a known Jxl column vector, and 
{P} is a column vector of discrete pressure values p(x); these matrices result 
from the evaluation of equation (4) . 

Assuming for the moment that the driving frequency, 03 , is not at (or very 
near) certain characteristic wave numbers of the fluid field enclosing the 
structure, equation (7) can be solved for {P}, thus 

{p> = [z(to)]{w> + [Lrhp 1 } (8) 


where [ Z (oj) ] = [L] ^[R] • 

Reference 13 has presented a method for arriving at equation (8) even in 
situations where CO is at or near one of the characteristic cavity resonance 
frequencies.* Briefly stated, the improved method consists of determining the 
unique surface pressure, p(x), that simultaneously satisfied the surface 
Helmholtz integral equation (4) and the interior Helmholtz integral reference 
14 . The interior Helmholtz integral is a relation similar to the form of 
equation (4) except it relates the fact that the fluid pressure for all points 
in the region of space occupied by the structure is zero. Enforcing this 
interior Helmholtz integral over a judiciously selected set of M interior points 
leads to a matrix analogous to equation (7) in the form 

= [R^iw} (9) 

Thus, equations (7) and (9) result in a set of (J+M) equations for the J unknowns 
{P}. Solving the overdetermined set of equations specified by equations (7) and 
(9) in a least square sense leads to an equation in the same form as equation (8) 
except that Z(co) is determined in a more involved manner. 

Next, by employing the principal of virtual work, the total surface 
pressures can be related to a set of consistent interaction nodal forces, {f}, 
thus 

{F} = [C I ]{P> 


* This is sometimes referred to as the cavity resonance problem. 
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( 10 ) 


For harmonic steady state problems, the continuous velocity, w(t), and dis- 
placement, u ( t) , amplitudes are related by w(t) = iwu(t). Making use of this 
relation in conjunction with the surface geometry relating normal components 
of motion into the Cartesian components employed in equation (1) results in 
the expression 


{w> = iu)[S]{u} 

Thus, combining equations (10), (8), and (11) leads to 

(F> = [T]{u} + [C I ][L]" 1 {P :L } 


(ID 


( 12 ) 


where [T] = iw[C*][Z][S] is typically a fully populated matrix that relates 
the interaction forces to the boundary displacement field. 

For steady state harmonic motion, all response qualities are proportional 
to exp(+iu)t). Thus, (u) = {u}°exp(itot) , { f} = {F)°exp(iO)t) and {F e ) = {F }° 
exp(itot) and the corresponding equation of motion for the structure (equation 
(1)) becomes 


(-<d 2 [M] + ia)[C] + [K])(U}° = -{F>° + (Fjj,) ° (13) 

where e"^ 1 " has been canceled out on both sides of the equation. Thus, substi- 
tuting equation (12) into equation (13) results in the relation 

[V]{U>° = {F } ° (14) 

A 

where [V] = -u 2 [M] + iw[C] + [K] + [T] and {F.}° - {fT 0 - [C 1 ] [Lj^tp 1 } 

A £« • 

It is to be noted that equation (13) contains matrices that are the size 
of the entire structure whereas the matrix size in equation (12) is only a size 
corresponding to the nodes in contact with the fluid. Thus, when substituting 
equation (12) into equation (13), allowances must be made in filling out the [T] 
and product matrix [C*] [L] -1 with zeros in the appropriate place to account for 
the matrix size mismatch. 

Formally, one may now state the solution to the interaction problem as 
finding the inverse of the highly populated [V] matrix. Thus, 

{U>° = [V] -1 (F }° (15) 

A 

Once {u}° is determined all other response quantities can be routinely computed. 
Substituting the solution {u}° into equation (6) and equation (11) and then 
equation (11) into equation (8) provides the total pressure, {p}, at the inter- 
face. Then substituting the surface pressure and surface velocity into the 
exterior form of the Helmholtz integral, reference 3, the pressure in any far 
field point in the media can easily be computed. Premultiplying the surface 
motion, {u}°, by the individual (unassembled) stiffness matrix for each element 
produces the individual structural nodal forces which in turn can be converted 
to element stresses. oqo 


For large size problems, the nonsymmetry and highly populated form of the 
complex [V] matrix makes its inversion somewhat of a problem when [V] is large. 
In some situations, [V] is ill-conditioned for certain frequency ranges due to 
the presence of large size [K] terms in the [V] expression in comparison to the 
rest of the terms comprising [V]. To get around these problems, an alternate 
modal analysis approach is sometimes taken, references 4 and 15. 

For the modal approach, let [(f)] be the NxM matrix of M undamped, in vacuo 
modes of the structural vibrations having N degrees of freedom, thus 

[*] = (16) 

th 

where (ij; (x)} is the m mode column vector which is normalized to the MxM unit 
identity matrix [I] such that 


[ 4 >] T [m] [ 4 .] = [i] (17) 

The modes [(f)] have the property that 

[<1>] T [K] [<j>] = [X] (18) 

where [X] is a MxM diagonal eigenvalue matrix whose non-zero elements are the 
squares of the natural frequencies (rad/sec) of the structure. The displace- 
ment field can be expressed in terms of the modes by the relation 

(u}° = [<J>]{Q}° (19) 

Next, upon substituting equation (19) into equation (14) and premultiplying the 
result by [(f)] T , one obtains 

[<t>] T [V][<j>]{Q}° = [<I>] T {f a } 0 (20) 

which can be rewritten in short notation as 

[V]{Q}° = {F g >° (21) 

where [V] = -co 2 [I] + [X] + [ 4> ] T ^iai [ C ] + [T]^[<J>] (22) 

and {F g )° E [<j>] T {F A }° (23) 

Generally, the MxM [V] matrix is complex, nonsymmetric and only under special 
situations is the [V] matrix fully diagonal (note only the first two contribu- 
tions to equation (22) are diagonal). When [V] is fully diagonal, its inver- 
sion is trivial, however, the general case must usually be considered where 
one is faced with the inversion of the [V] matrix in order to solve the system 
of equations defined by equation (21). Formally, then, the solution to the 
fluid structure interaction problem can be written as 

(q}° = [V] - 1 {F g }° 
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(24) 


where we have traded having to invert a NxN [V] matrix in the direct approach 
for having to invert MxM [V] matrix in the modal approach. Strictly speaking, 
there is one mode shape for each degree of freedom, consequently, if M is set 
equal to N, one is right back where one started in being faced with the inver- 
sion of a NxN complex matrix. However, one can usually judiciously relate the 
important modes of vibrations based on certain symmetries of loading or based 
on the customary omission of the higher modes of vibration. After the selection 
process, one is usually left with a [V] matrix that is substantially smaller in 
size than the original [V] matrix encountered in the direct approach. 


The previous development is for a general shaped submerged structure. 

In the special case where the body has an axis of revolution, it is possible 
to use a Fourier series decomposition in the angular variable of a cylindrical 
coordinate system centered about the axis of revolution. Thus, one can describe 
an arbitrary pressure (or velocity) distribution through the relations 


p(x,6) 



p (x)Cosn0 + 
n 



p(x)Sinn0 


w(x,0) 


CO 



w (x)Cosn0 + 
n 


n=l 


w (x)Sinn0 
n 


(25) 


Applying such an expansion to the development just presented for the general 
three-dimensional case results in problem formulation analogous to equation (14) 
or to equation (21). The main difference is that the coefficient matrix [V(n)] 
in equation (15) (or corresponding [V(n)] matrix in equation (21)) is now a 
function of the wave number n corresponding to the expansions in equation (25). 
Consequently, solving the full three-dimensional problem is equivalent to solving 
a sequence of n = 1, 2, . . N smaller sets of linear equations (i.e equations 
(14) or (21)). The phrase "smaller sets of linear equations" is used since the 
elimination of the third spacial dimension (through the introduction of the 
Fourier expansion) substantially reduces the size of the coefficient matrix [V]. 
From a computational point of view, it is usually more efficient to solve, say, 
six (H=6) two-dimensional size problems than one large three-dimensional one. 

The details of setting up the actual [V(n)] array, for cubic polynomial dis- 
placement fields, is presented in more detail in reference 4. Consequently, 
it will not be repeated here. 


The interface of this solution technique with the NASTRAN computer 
program can be made in one of two ways where the selected approach depends on 
whether the direct or modal solution technique is used to solve the problem. 

In the case of the modal formulation, a computer program called FIST (Fluids 
Interacting with Structures) has been written which directly accepts, as input 
from NASTRAN, the structural mode shapes for the in vacuo normal modes or mass 
and stiffness matrices. With the addition of a few simple alter cards in the 
NASTRAN run stream, the modes are written on tape from NASTRAN using the module 
OUTPUT 2. For example, in rigid format 3, adding the cards 
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ALTER 84 

OUTPUT 2 KAA,MAA,,,//C,N,1/C,N,0 $ 

ALTER 94 

OUTPUT 2 PHIA,MI,,,//C,N,1/C,N,0 $ 

ALTER 96 

OUTPUT2 PHIG,, , ,//C,N,l/C,N,0 $ 

ENDALTER 

in the EXECUTIVE control deck is all that is needed to have the information 
needed by FIST to solve the flu id /structure interaction problem. The current 
version of FIST is currently limited to axisymmetric structures subject to non- 
axisymmetric loadings. 

In the case of the direct formulation, FIST has an option which permits 
one to solve the system of linear complex simultaneous equations (14). However, 
the current version is limited to a 40 x 40 nonsymmetric , fully populated 
matrix. In situations where the [V] matrix is large, the NASTRAN mathematical 
solution routines can be utilized to solve the problem at hand. The frequency 
response rigid format number 8 solves the following problem 

j^[M]0) 2 + iU)[B] + [K]J{X> = {P} (26) 

for the displacement amplitude {x}, where the multiplying matrix on {x} can 

be, in general, nonsymmetric and complex. The complex loading vector, {p), 

can be introduced into NASTRAN through the RL0AD1 bulk data card and the [M] , 

[ B ] , and [K] arrays need not be computed by NASTRAN but rather can be directly 

inserted, element by element, through a DMIG card. Since no structure is given 

to NASTRAN directly, all the [M] , [B] and [K] matrices are zero, except for the 

direct input components (denoted by [M 2 ;,], [B^J, and [K* ] in the NASTRAN 

da dd dd 

theoretical manual, pg. 9.3-7). The [M] and [B] arrays are zero by virtue of 
not defining them in any way. Thus, there remains the [K* ] array which is 
read in (in complex form) via the DMIG card. NASTRAN willproceed in the 
usual manner for the direct frequency response solution and compute the (x) 
solution which, of course, corresponds to the desired result (u}° (i.e. equation 
(15)). 


By employing alter instructions, one who is familiar with DMAP can most 
likely perform the desired operations in a more direct fashion. The advantage 
of the dummy stiffness application method is that is can all be done within 
the current fixed format of the program. 

As an illustration of the solution technique, consider the situation 
where a steel submerged thin wall spherical shell is harmonically driven by a 
point concentrated force. The shell has a radius of 2.54 cm, wall thickness 
of 0.127 cm and a nondimens ional driving frequency of KA = 0.4 (K 5 u)/c where 
U) is the driving frequency in rad/sec and c is dilatational wave speed in the 
water). The exact solution to this problem is shown in figure 1 by the solid 
line (reference 12) and is processed using the first 50 axisymmetric modes of 
a thin sphere. The fine dashed line corresponds to a FIST solution to the 
problem (employing a cubic distribution pressure variation) which uses the 
direct solution approach (equation (15)) with the sphere subdivided into eight 
segments. The coarse dashed line corresponds to a modal solution (equation (24)) 
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that employs three mode shapes and an eighteen segment subdivision. Both FIST 
solutions employ too coarse of a mesh to predict an accurate response near the 
point load. Agreement away from the point load is seen to be very good. 


STEADY STATE-COUPLED FLUID FIELD 


Again, the case of a totally submerged elastic structure is treated, 
except that in this case the structure as well as the fluid is modeled with 
finite elements. The first approach, of the two presented in this subsection, 
is the situation where displacement (or "mock") fluid elements is employed in 
the problem formulation, references 17, 18. In this situation, we start with 
the three-dimensional elasticity equations 


8V 
P ItT 


!!±i 

9x. 

J 


i = 1,2,3 


(27) 


with O . . 
13 



3x t 

k 


+ y 



+ 




is the Kronecker delta) 


(27a) 


We start by letting the Lame constant y-H); the Lame constant A+k; and noting 
that the reduced equation (27a) now implies that 


°11 a 22 °33 


and a 


12 


= a 


31 


= a 


32 


= a 


21 


°13 ' a 23 ‘ °' 


(28) 


Finally, upon defining P = -a.^ = “022 = -0 33 ’ ^ i s seen that the reduced 
equations (27) and (27a) represent tne same^field equations as those defined 
by equations (2, 2a). 

Consequently, any solid elements in NASTRAN that are built from working 
with equations (27, 27a) can be converted into mock fluid elements by appro- 
priately refining the constants in the elasticity stress-strain law 

{a} = [G]{e} 


For three-dimensional type elements like brick and ring elements, the array of 
elastic constants for an isotropic material can be written as 


U+2y) 

A 

A 

0 

0 

0 

A 

(A+2y) 

A 

0 

0 

0 

X 

A 

(A+2y) 

0 

0 

0 

0 

0 

0 

y 

0 

0 

0 

0 

0 

0 

y 

0 

0 

0 

0 

0 

0 

y 


(29) 
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Similarly, for two-dimensional solid elements, like membrane elements, the 

array of elastic constants (corresponding to {a , a , a } can be written as 

xx yy xy 

G 11 G 12 G 13 

G 21 G 22 G 23 

G 31 G 32 G 33 

For the NASTRAN program installation of the mock elements for 2-D membrane 
solid elements is achieved through a MAT2 card wherein 

G = G 12 = G 21 = G 22 = k S fluid bulk modulus 
and the remaining G^ = G^ = G^^ = G^ 2 = G^ = 0.0. 

For employing three-dimensional mock elements, the situation would be 
straightforward if the [G] matrix were allowed to be input in a general form 
like the two-dimensional case; or, if the stress-strain matrix were written in 
terms of Lame constants rather than in terms of the more common modulus of 
elasticity, E, and Poisson's ratio, V. For the latter case, one need only set 
y = 0 and X = k and the desired mock element could be formed. In actuality, 
the NASTRAN [G] array accepts input in the fori# of E and V and internal to the 
program, the elements of the [G] array are defined as follows 

G 11 = G 22 = G 33 = *<1-V)/[(1-2V)U4V)] 

G 12 ' G 21 " G 13 ' G 31 ' G 23 ' G 32 ' 1 (31) 

G 44 ' G 55 * °66 ‘ • 5E/(1+V) 

The problem can be resolved by rewriting a small portion of the Fortran 
coding that fills out the material constants array in the desired form. 

An alternate procedure is to adjust the values of E, V (or G) on a MAT1 
card so that X = k and y = 0. Setting G = 0 (note that y = G) and solving for 
the E required to be consistent with the proper fluid bulk modulus, k, will not 
work because the NASTRAN coding tries to form V by dividing G (see footnote) 
and dividing by zero will not be handled properly by the computer. 

* The shear modulus, G, can be given in place of E (or V) in which case the E 
or V is computed internally from E = 2G(l+v) (or E 1 ). 

V “ 2G 
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One way to resolve the problem is to relax the strict equalities for 
the mock element that A = k and p = 0 but rather enforce them only approximately 
such that A~k and y~0 (i.e. so long as A»>]j). By setting 


and 


V = 0.49999 


k(l-2v) (1+v) 


k/16,664.44 


Ola) 


on the MAT1 cards, NASTRAN will internally generate a set of . constants that 
will adequately represent (but not be exactly equal to) the desired exact values. 
As an illustration, consider water that has a bulk modulus of k = 316,000 psi 
(2.18xl0 9 N/m 2 ). Thus, for p = 0, typical proper values of the [G] array would 
correspond to 


G n = 316,000, G 12 = 316,000, and G ^ = 0.0. 

Applying the suggested approximate approach via equations (31a) , NASTRAN invokes 
equations (31); thus would internally compute the typical [G] array elements as 

G = 316,050 G 12 = 316,037 G ^ = 6.32 

which should be sufficiently close to produce fluid response results of the 
same degree of accuracy had the exact [G] entries been used. 

Pressure distribution information is obtained by examining the stress 
output from NASTRAN. The fluid pressure is obtained by reversing the sign of 
the normal stress output (since all the normal components are equal for mock 
elements, the user can select any normal component). 

The boundary between the fluid and solid is handled by only forcing the 
normal component of fluid displacement to be equal to the normal displacement 
of the interfacing solid. This can be easily done through the introduction 
of a double node in conjunction with a MPC constraint. The boundary at 
infinity is handled analogous to the approach used in solids for earthquake 
problems (reference 19) and later in fluid applications (references 17, 18). 

This is accomplished by placing the fluid boundary not at infinity, but at a 
finite distance that is far enough so that interaction waves radiating from 
the submerged structure will satisfy (or nearly satisfy) the boundary condition 

P = Pcu R (32) 

where p is the fluid mass density, c is the fluid sound speed and u is the 
velocity normal to the outer boundary. This condition is true for plane waves 
and asymptomatically true for cylindrical and spherical waves. The finite 
element form of equation (32) is given by 

{F fe } = [C b Hu> (33) 
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where [C^] is a diagonal matrix with zero diagonal values for non-outer 
boundarypoints and a value of pcAA. for normal outer boundary degrees- of- 
freedom where AA. is the pressure- to-f or ce conversion term and corresponds to 
a segment area at the boundary node. In NASTRAN, the boundary dash pots are 
applied with CDAMP1 cards. 

A rough guideline is needed for determining how far the fluid boundary 
should be placed in order for the plane wave approximation, equation (32), to 
be valid. A steady state solution can be constructed from some distribution 
of point sources around the structure fluid boundary. For a 3-D problem, a 
single source will approach (within 98.6%) a plane wave pressure velocity 
relationship (like equation (32)) after moving one wave length away from the 
source. # The percentage quoted refers to the fact that complex impedance 
(z = p/u is .986pc. Moving away 11/2 wave lengths, this percentage becomes 
99.4%. Fhus, we are suggesting that the single source decay information can 
be used to judge the distance to place the absorbing boundary. For 2-D 
problems, a line source emitting cylindrical waves will approach a plane wave 
boundary condition to within 99.2% for one wave length away and to within 99.7% 
for 11/2 wave lengths away. Thus, it is suggested to place the boundary a 
distance D away from the structure where 

D = aA (34) 

and a is a proportionality constant (e.g. 1.5 for a 99.4% correct plane wave 
assumption} related to the degree of the plane wave boundary condition assump- 
tion, and A is the wave length of the steady state driving frequency in water 
(i. e. X = 2ttc/u)) . 

Next, one must consider the size elements to use so that the elements of 
the mesh do not artificially "ring 11 at their natural frequencies. To avoid 
ringing, there exists a minimum element length, A L , that is related to some 
fraction, 6, of the wave lenght of the driving frequency in the fluid, thus 

AL = 3A (35) 

The value of 8 will depend on the type of elements being used. For example, 

if one employs CQDMEM elements of the NASTRAN program, 8 ~ 1/6 to avoid mesh 

ringing. Modeling the region from the fluid structure interface out to the 

mathematical cut in the fluid boundary would result in n elements of length 

A , thus 
L 



Substituting equations (34) and (35) into equation (36) results in the expression 

a 

n “ 8 

which is independent of the driving frequency 0). Thus, employing typical values 
of a = 1.5 and 8 - 1/6 into equation (36) shows, for example, that regardless 
of the driving frequency magnitude, it is possible to model the fluid field with 
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as few as 9 elements in the direction normal to the surface. In cases where 
the surface structure elements are coarse, more elements would be needed to 
blend in the fine surface elements into the coarser field elements. 

For radiation type problems (i.e. a structure vibrating and transmitting 
outgoing waves) , one can handle the infinite boundary problem by placing the 
dash pots around the outer boundary. For scattering type problems where a 
submerged structure is subject to an incident wave input, the handling of the 
infinite boundary problem is more complicated in that both the incident wave 
driving force and the boundary dash pots both appear on the boundary. The 
manner in which this class problem is handled is discussed in reference 17 and 
will not be repeated here. 


As a demonstration problem consider the two-dimensional problem of an 
infinitely long cylindrical inclusion imbedded in an acoustical fluid medium. 
The inclusion is subjected to an incident harmonic plane wave (this is classi- 
fied as a scattering problem) . A sketch of the full model and corresponding 
finite element sketch is illustrated in figure 2. Constant strain CQDMEM 
elements are used to construct the model for both the solid and mock fluid 
elements. Two solutions (for two different inclusion types) are presented in 
the form of a comparison between an exact and corresponding finite element 
response. In either case, the exact solution (references 20, 21) is repre- 
sented by the solid curve and the dots are the corresponding finite element 
solution. The solution response is given, in non-dimensional form, as the 
ratio between the total pressure to the incident free field amplitude. The 
parameters in the upper right corner of the figure denote a set of non-dimen- 
sional parameters that characterize the physical parameters of the problem and 
have the corresponding definitions 


T /- x __ i j ^ ^inclusion radius 

KA = non-dimensional driving frequency = -tt — r*; ; — 

fluid wave speed 

r - radial coordinate 
inclusion radius 

5 _ fluid dilatational wave speed 

inclusion dilatational wave speed 

~ _ inclusion dilatational wave speed 

inclusion shear wave speed 

- _ fluid mass density 

inclusion mass density 


The response shown in figure 3 corresponds to a vacuous inclusion and the 
response in figure 4 to an elastic aluminum inclusion. Except for the 0° and 
180° (back and front) data points on the aluminum cylinder, the response results 
agreement was good. Response comparison for other radii (both closer and 
further away from the results presented) gave equally good results. The mesh 
size used was pushing the limit regarding the size needed to avoid ringing. 

It is felt that a finer mesh would have improved the results in the 0° or 180° 
data points for the aluminum solutions. 
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Reference 22 presents another completely different approach to solving 
steady state problems with NASTRAN. In this approach the submerged structure 
is surrounded with a sphere shaped region of finite elements. The proper 
boundary condition for handling the infinite fluid domain beyond the bounding 
sphere surface is treated with an eigenvalue expansion approach. The details 

are presented in reference 22 (paper in this colloquium), therefore will not be 
repeated here. 


TRANSIENT-UNCOUPLED FLUID FIELD 

This class of problem treats a totally submerged structure subject to 
dynamic loading usually in the form of an incident pressure wave. Reference 23 
presents an application of the NASTRAN program towards solving transient fluid/ 
structure interaction problems with the DAA (doubly asymptotic approximation) 
method. This method involves imposing an analytical decoupling relation (in 
differential equation form) describing the relationship between the pressure 
at the interface and the corresponding interface motion. The decoupling 
approximation eliminates the need for modeling the fluid field with finite 
elements. A brief outline of the method is presented in reference 23. A more 
detailed discussion of the impl mentation into NASTRAN is presented as a paper 
in this colloquium (reference 24), therefore the reader is referred to that 
paper for more details. 


TRANSIENT-COUPLED FLUID FIELD 


The class of problem considered here is the same as the previous transient 
category except in this case the fluid is modeled as part of the finite element 
network. The first NASTRAN application in this category considers the case 
where pressure type fluid elements are used to model the fluid field. These 
types of elements are different from the displacement type elements discussed 
earlier in that there is only one degree of freedom per node (namely pressure); 
this is in constrast* to one, two or three degrees of freedom per node for 
mock fluid elements which have displacements as the basic unknowns. The imple- 
mentation of this method into NASTRAN requires one to dummy the construction 
of the stiffness and mass matrices of a conventional displacement type finite 
element so that only one displacement component is active (the remaining ones 
are zero); and further, the remaining nonzero component plays the role of 
pressure. The proper units are handled through redefining the elements of the 
t G ] stress-strain matrix. This approach to solving fluid/structure interaction 
problems with NASTRAN was first introduced in the 4th NASTRAN User’s Colloquium 
(reference 23) . The complete details of the implementation of the method is 
presented in the current 5th NASTRAN User’s Colloquium (reference 24). 

The second application of NASTRAN (employing the coupled fluid field 
approach) is that of using the mock fluid elements. These displacement type 
elements were already discussed in the previous secion on harmonic analysis. 

The method of implementing them via the stress-strain matrix [G] is done in 


* The actual number depends on whether one is solving a one, two, or three- 
dimensional problem. 
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exactly the same manner as for the transient type problems as well. Transient 
problems using mock elements are solved using rigid format number 9. The 
infinite fluid boundary problem can be handled by temporal truncation. This 
is the most straightforward approach and is readily adaptable to both the 
pressure or mock element type transient solutions. In construction of the 
finite element mesh, one models the fluid surrounding the solid cutout to, 
say, 2 structure lengths. For scattering or radiation type problems, one can 
take advantage of the fact that the continuous equations are hyperbolic in 
nature. Thus for, say, a radiation problem, the solution will be such that the 
response in front of a radiating wave is zero, thus the problem does not know 
a boundary to the mesh even exists until the radiating wave actually gets there 
(due to the discretization of the problem, the governing equations do not 
exactly behave like hyperbolic equations, but the idea of traveling waves are 
still roughly approximated by the governing discretized differential equations). 
Thus, the solution to the problem can be obtained in the same manner as a 
finite boundary case, except that the solution response must be truncated at 
the time when the radiated wave reaches the mesh boundary. 

Scattering problems can be treated in a similar manner. The free-field 
incident wave solution starts the problem, i.e. the initial conditions for the 
problem solution are obtained by setting the response field, behind the inci- 
dent wave, equal to the free-field solution. The equations of motion are 
integrated in time in the usual manner but must be truncated when scattered 
waves off the structure reach the finite fluid mesh boundary. 

One-dimensional wave propagation examples using mock elements in NASTRAN 
are presented in reference 17*. A solution employing mock elements for a three- 
dimensional problem is shown in figure 5. The problem corresponds to a plane 
step wave traveling through a solid elastic homogeneous medium and interacting 
with a spherical cavity filled with fluid. At the time these runs were made, 
a level of NASTRAN capable of solving axis of revolution solid strcutures 
subject to nonsymmetric loading was not yet available, consequently, a pro- 
gram (reference 25) other than NASTRAN was used to generate the result shown 
in figure 6. The program uses a harmonic decomposition in the angular coordinate 
of a cylindrical coordinate system to reduce the full three-dimensional problem 
to that of solving a set of smaller two-dimensional ones (with the harmonic 
wave number as a parameter in each two-dimensional subproblem) . The full three- 
dimensional response is obtained by superposition over the angular harmonics 
(usually 5 terms are adequate). The response shown in figure 6 is the pressure 
at the center of the fluid sphere for both the exact solution (reference 26) and 
the corresponding finite element solution. For reference, the pressure in the 
free field (negative of the average normal stress in the solid) is also shown. 
Transient finite element solutions of this type tend to ring about the true 
solution. The frequency of the ringing is associated with the highest natural 
frequency of the mesh. Reference 27 discusses this point in detail and provides 
a digital filtering technique for eliminating some of the ringing problem. 


* In reference 17, there are several sign errors that should be pointed out to 
avoid confusion with the development presented here; namely, on page 74 of 
reference 17, replace k with -k in equation (3), replace -k with +k on the 
third line from bottom and , finally, replace -k with +k at the bottom of page 75. 
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For completeness, a brief discussion of the direct application of the 
NASTRAN computer program for solving fluid/structure interaction problems is 
presented next. The NASTRAN program has built-in pressure type elements (called 
CFLUID1 elements) and are described through RINGFL , PRESPT and FREEPT type fluid 
nodes. The elements are designed to operate in contained tanks that may have 
either rigid or elastic walls. These special elements have the following 
restrictions: 

1. The user may not apply loads, constraints, sequencing or omitted 
coordinate directly on the fluid nodes involved. Instead, the user supplies 
information related to the boundaries and NASTRAN internally generates the 
correct constraints, sequencing, and matrix terms. 

2. The input data to NASTRAN may include all of the existing operations 
except the axisyrametric structural element data (e.g., ax i symmetric shell 
elements cannot be used) . 

3. The fluid must lie within the walls of an open or closed tank. 

4. The first 6 rigid formats of NASTRAN may not be used in conjunction 
with these elements. NASTRAN assumes the walls of the container are rigid for 
these first 6 rigid formats but allows elasticity for the remaining 6 (fortunately 
direct frequency and direct transient response are included in the remaining 6) . 

5. No means are provided for the direct input of applied loads on the 
fluid. Loading must come through the motion of the walls. 

The list of constraints that are placed on the usage of these elements 
rather severely limits the range of application, particularly in the case where 
unbounded fluid regions are of interest. Even within these constraints (refer- 
ence 28), however, some rather interesting applications to acoustic noise problems 
associated with automobiles have been found. 

Another unique feature that the NASTRAN program has is the ability to 
treat the free surface problem and include gravity terms into the fluid equations 
of motion. 


CONCLUSIONS 

The latest version of the NASTRAN computer program, as of this writing, 

does not handle a very large class of fluid/structure interaction problems via 

the direct rigid format application of the program. This paper presents a 

variety of nonstandard usage of the program to broaden the scope of problem 

application in the area of fluid/structure interaction. The implementation 

of the techniques presented here varies from one extreme of requiring the user 

to only make slight modifications to the standard problem input of NASTRAN-to- 

another extreme of requiring a substantially sized auxiliary support computer 

program to handle pre- and/or post-processing of the input and output data. 

The implementation of these methods depends, to some degree, upon the ingenuity 

of the user. Hopefully, future versions of NASTRAN will have more automatic 

procedures for solving problems of the type addressed in this paper. 

* 
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SURFACE PRESSURE = P*(4irA 2 )/F 


EXACT SOLUTION 











1 . 


EXACT 

(7) FINITE 
W ELEMENT 

KA = 2.122 

R = 2.125 


Figure 3. - Total (Incident & Scattered) Steady State 
Pressure Response of a Cylindrical Void in Water 


EXACT 

(7) FINITE 
ELEMENT 

KA = 2.122 

R = 2.125 

5 = 0.2445 

C = 1.980 

p = 0.380 


Figure 4. - Total (Incident & Scattered) Steady State Pressure 
Response of a Solid Elastic Aluminum Cylindrical Inclusion In Water 
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CENTER PRESSURE r (STATIC PRESSURE 
RESPONSE) 


AXIS OF REVOLUTION 



2.q 


O FREE FIELD INPUT 
• EXACT 

A FINITE ELEMENT 



p g = SOLID DENSITY 

p = FLUID DENSITY 

V = SOLID POISSON'S RATIO 

C = SOLID DILATATIONAL WAVE 


A = CAVITY RADIUS 


Figure 6. - Pressure Transient Response at Fluid Cavity Center 
(p / p. = 1; C /C = 1; V = 1/4) 

SI SI 
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